C***************************************************************
C      Subroutine L2LC : www.boundary-element-method.com                       
C***************************************************************
C 
C  Version 1.1 (Sept 2015) Inclusion of revised GEOM2D file.
C   Geometry verification now in VG2LC, and in-built LOG and  
C   SQRT functions used.
C  [Version 1 2001]
C  Stephen Kirkup, 
C  http://www.researchgate.net/profile/Stephen_Kirkup
C  School of Engineering 
C  University of Central Lancashire 
C
C  This open source code can be found at
C   www.boundary-element-method.com/fortran/L2LC.FOR
C  The use guide for the L2LC code can be found at
C   www.boundary-element-method.com/fortran/L2LC_FOR.pdf
C  A test harness for this code can be found at
C   www.boundary-element-method.com/fortran/L2LC_TESTS.FOR
C
C  Issued under the GNU General Public License 2007, see gpl.txt
C
C***************************************************************
C This subroutine computes the discrete form of the 2-dimensional
C Laplace integral operators L, M, Mt, and N over a straight line
C element. The subroutine is useful in integral equation methods for
C the solution of 2-dimensional Laplace problems.
C
C The following diagram shows how the subroutine is to be used. A main
C program is required along with the definition of  LOG and
C SQRT. The package is the contents of this file.
C
C                                   ....................................
C                                   :                                  :
C                                   :                                  :
C      ----------------------       :     --------------------------   :
C      |                    |       :     |                        |   :
C      |   MAIN PROGRAM     |------->-----|           L2LC         |   :
C      |                    |       :     |                        |   :
C      ----------------------       :     --------------------------   :
C                                   :                 |                :
C                                   :                 >                :
C                                   :                 |                :
C                                   :      ------------------------    :
C          Package ---------------->:      | Subordinate routines |    :
C                                   :      ------------------------    :
C                                   :                                  :
C                                   :..................................:
C
C The subroutine has the form:
C      SUBROUTINE L2LC(P,VECP,QA,QB,LPONEL,
C     * MAXNQ,NQ,AQ,WQ,
C     * LVALID,EGEOM,EQRULE,LFAIL,
C     * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)
C
C The parameters to the subroutine
C ================================
C Point (input)
C real P(2): The Cartesian coordinates of the point p.
C real VECP(2): The Cartesian components of the unit normal at p. 
C  Required in the computation of DISMT and DISN. The squares of the
C  components must sum to one.

C Geometry of element (input)
C real QA(2): The Cartesian coordinates of the first edge of the element.
C real QB(2): The Cartesian coordinates of the second edge of the element.
C logical LPONEL: If the point P(2) lies on QA-QB then LPONEL must be set
C  .TRUE., otherwise LPONEL must be set .FALSE..

C Quadrature rule (input)
C integer MAXNQ: The limit on the size of the quadrature rule. The value
C  should not be changed between calls of L2LC. MAXNQ>=1.
C integer NQ: The actual number of quadrature rule points. 1=<NQ<=MAXNQ.
C real AQ(MAXNQ): The quadrature rule abscissae. The values must lie in
C  the domain [0,1] and be in ascending order.
C real WQ(MAXNQ): The quadrature rule weights which correspond to the
C  quadrature points in AQ. The components of WQ must sum to one.

C Validation and control parameters (input)
C logical LVALID: A switch to enable choice of checking of subroutine
C  parameters (.TRUE.) or not (.FALSE.).
C real EGEOM: The maximum absolute error in the parameters that
C  describe the geometry. Value is of importance only when LVALID=.TRUE..
C real EQRULE: The maximum absolute error in the components of the
C  quadrature rule data. Value is of importance only when LVALID=.TRUE..

C Validation parameter (output)
C logical LFAIL: Value is only important if LVALID=.TRUE.. If 
C  LFAIL=.FALSE. then the input data has been found to be satisfactory. 
C  If LFAIL=.TRUE. then the  input data has been found to be 
C  unsatisfactory. The subroutine would have been aborted. The output
C  parameters DISL, DISM, DISMT and DISN will all be zero. A 
C  diagnosis will be given in the file L2LC.ERR.
C Choice of discrete forms required (input)
C logical LL: If discrete form of Lk operator is required then set
C  .TRUE., otherwise set .FALSE..
C logical LM: If discrete form of Mk operator is required then set
C  .TRUE., otherwise set .FALSE..
C logical LMT: If discrete form of Mkt operator is required then set
C  .TRUE., otherwise set .FALSE..
C logical LN: If discrete form of Nk operator is required then set
C  .TRUE., otherwise set .FALSE..

C Discrete Laplace integral operators (output)
C complex DISL: The discrete Lk integral operator.
C complex DISM: The discrete Mk integral operator.
C complex DISMT: The discrete Mkt integral operator.
C complex DISN: The discrete Nk integral operator.

C External modules to be supplied
C ===============================
C real function SQRT(X): real X : Must return the square root of X.
C real function LOG(X): real X : Must return the natural logarithm
C  of X.
 
C Notes on the validation parameters
C ----------------------------------
C  (1) If LVALID=.TRUE. then a file named L2LC.ERR should be open when
C   subroutine L2LC is entered. Use a command such as 
C         OPEN(UNIT=10,FILE='L2LC.ERR',STATUS='UNKNOWN')
C  at the beginning of the calling program and the corresponding command
C         CLOSE(10)
C  at the end of the calling program.
C  This file accumulates the error messages and warnings from the 
C   subroutine. If this file is not open then an error message is 
C   output to standard output and LFAIL=.TRUE.
C (2) If LVALID=.TRUE. then EGEOM must be less than 1% of the maximum
C  absolute coordinate of QA and QB.
C (3) If LVALID=.TRUE. then EQRULE must be less than 10%/NQ.

C Notes on the geometric parameters
C ---------------------------------
C  (1) P, QA and QB must be distinct points (with respect to EGEOM).
C  (2) If LPONEL=.TRUE. then P must lie on element QA-QB. If 
C   LPONEL=.FALSE. then P must not lie on QA-QB.
C  (3) The normal to the element is defined to be
C             [-QB(2)+QA(2), QB(1)-QA(1)] normalised.
C   Hence the normal is to the right on the line QA-QB.
C  (4) If LPONEL=.TRUE. (and either LMT=.TRUE. or LN=.TRUE.) then 
C   VECP must be [-QB(2)+QA(2), QB(1)-QA(1)] normalised.

C Notes on the quadrature rule parameters
C ---------------------------------------
C  (1) The quadrature rule is assumed to compute linear functions
C   exactly (with respect to EQRULE).
C  (2) If LPONEL=.TRUE. then the first derivative of the Lk and Nk
C   functions on the element are discontinuous at P. The input
C   quadrature rule should take account of this by using a
C   composite rule that divides at P in these cases.

C External modules provided in the file GEOM2D.FOR
C ================================================
C real function SIZE2: Returns the modulus of a 2-vector.
C real function SSIZE2: Returns the square of the modulus of a 2-vector.
C real function DOT2: Returns the dot product of two 2-vectors.
C Subroutine VEC2: Gives the output in VEC for the subtraction.

C External module provided in the file VG2LC.FOR
C ================================================
C Subroutine VG2LC: Verifies the geometrical input to L2LC.

      SUBROUTINE L2LC(P,VECP,QA,QB,LPONEL,
     * MAXNQ,NQ,AQ,WQ,
     * LVALID,EGEOM,EQRULE,LFAIL,
     * LL,LM,LMT,LN,DISL,DISM,DISMT,DISN)


C Point
      REAL*8     P(2)
      REAL*8     VECP(2)

C Geometry of element
      REAL*8     QA(2)
      REAL*8     QB(2)
      LOGICAL    LPONEL

C Quadrature rule
      INTEGER    MAXNQ
      INTEGER    NQ
      REAL*8     AQ(MAXNQ)
      REAL*8     WQ(MAXNQ)

C Validation and control parameters
      LOGICAL    LVALID
      REAL*8     EGEOM
      REAL*8     EQRULE
      LOGICAL    LFAIL

C Choice of discrete forms required
      LOGICAL    LL
      LOGICAL    LM
      LOGICAL    LMT
      LOGICAL    LN

C Discrete Laplace integral operators
      REAL*8     DISL
      REAL*8     DISM
      REAL*8     DISMT
      REAL*8     DISN


C External functions to be supplied
      LOGICAL    VG2LC
      LOGICAL    VQUAD

C External functions provided
      REAL*8     SIZE2
      REAL*8     SSIZE2
      REAL*8     DOT2

C Variables that remain constant throughout subroutine, once initialised
C  Constants
      REAL*8     ZERO,ONE,TWO,PI,TWOPI,OO2PI
C  L2LC.ERR open status and unit number
      LOGICAL    LOPEN
      INTEGER    IU
C  Useful geometric values
      REAL*8     QBMA(2),PMA(2),PMB(2),NORMQ(2)
      REAL*8     QLEN,PQALEN,PQBLEN,DNPNQ
C  Control values relating to required solutions     
      LOGICAL    LLORN,LMORN,LMTORN,LMSORN,LMORMT,LMSV,LMTSV

C Variables used in the validation of the input data
C  Error variables
      LOGICAL    LERROR

C Variables used in the operation of the quadrature rule
C  Variables for the accumulation of integrals
      REAL*8     SUML,SUMM,SUMMT,SUMN
C  Index of quadrature point
      INTEGER*4  IQ
C  Geometric variables
      REAL*8     Q(2),RR(2)
      REAL*8     SR,R,RNQ,RNP,RNPRNQ,RNPNQ
C  Green's functions for the Laplace operators and its r-derivatives
      REAL*8     G,GR,GRR
C  Other variables
      REAL*8     WGR


C INITIALISATION

C Set constants
      ZERO=0.0D0
      ONE=1.0D0
      TWO=2.0D0
      PI=3.1415926535897932
      TWOPI=TWO*PI
      OO2PI=ONE/TWOPI

C Initialise useful geometrical information
      CALL VEC2(QA,QB,QBMA)
      CALL VEC2(QA,P,PMA)
      CALL VEC2(QB,P,PMB)
      QLEN=SIZE2(QBMA)
      PQALEN=SIZE2(PMA)
      PQBLEN=SIZE2(PMB)


C Set LFAIL,LERROR to FALSE
      LFAIL=.FALSE.
      LERROR=.FALSE.

C If LVALID is false then bypass the validation process
      IF (.NOT.LVALID) GOTO 900


C VALIDATION OF INPUT

C Check that the file L2LC.ERR has been opened, if not then exit
      INQUIRE(FILE='L2LC.ERR',OPENED=LOPEN)
      IF (.NOT.LOPEN) THEN
        WRITE(*,*) 'ERROR(L2LC) - File for error messages "L2LC.ERR"'
        WRITE(*,*) ' is not open'
        WRITE(*,*) 'Execution of L2LC is aborted'
        LERROR=.TRUE.
        GOTO 998
      END IF

C Find out the unit number of the file L2LC.ERR
      INQUIRE(FILE='L2LC.ERR',NUMBER=IU)

C Set LERROR to .FALSE.
      LERROR=.FALSE.


      IF (EQRULE.LE.ZERO) THEN
        WRITE(IU,*) 'ERROR(L2LC) - Parameter EQRULE must be positive'
        LERROR=.TRUE.
      END IF
C  If LERROR then exit L2LC
      IF (LERROR) THEN
        WRITE(IU,*) 'Execution of L2LC is aborted'
        GOTO 999
      END IF

C Verify geometrical information
      LERROR=.NOT.VG2LC(P,VECP,QA,QB,LPONEL,EGEOM,IU)

C  If LERROR then exit L2LC
      IF (LERROR) THEN
        WRITE(IU,*) 'Execution of L2LC is aborted'
        GOTO 998
      END IF

      
C Verify quadrature rule
      LERROR=.NOT.VQUAD(MAXNQ, NQ, WQ, AQ, -1.0D0, EQRULE, IU)

C  If LERROR then exit L2LC
      IF (LERROR) THEN
        WRITE(IU,*) 'Execution of L2LC is aborted'
        GOTO 998
      END IF

900   CONTINUE


C FURTHER INITIALISATION

C Set the values  of the discrete operators to zero.
      DISL=ZERO
      DISM=ZERO
      DISMT=ZERO
      DISN=ZERO

C Exit if no solutions are required
      IF (.NOT.(LL.OR.LM.OR.LMT.OR.LN)) GOTO 998

C Initialise useful control information
      LLORN=LL.OR.LN
      LMORN=LM.OR.LN
      LMSORN=LM.OR.LMT.OR.LN
      LMTORN=LMT.OR.LN
      LMORMT=LM.OR.LMT
      IF (LPONEL) THEN
        LMSV=LM
        LMTSV=LMT
        LM=.FALSE.
        LMT=.FALSE.
      END IF

C Initialise useful geometric information
      IF (LMORN) THEN
        NORMQ(1)=-QBMA(2)/QLEN
        NORMQ(2)=QBMA(1)/QLEN
      END IF
      IF (LN) DNPNQ=DOT2(VECP,NORMQ)


C Initialise SUML-SUMN, the variables that accumulate 
C  the integrals Lk,Mk,Mkt,Nk (divided by QLEN)
      SUML=ZERO
      SUMM=ZERO
      SUMMT=ZERO
      SUMN=ZERO
      IF (LPONEL) THEN
        IF (LLORN) SUML=(ONE-(PQALEN*LOG(PQALEN)
     *   +PQBLEN*LOG(PQBLEN))/QLEN)/TWOPI
        IF (LN) SUMN=-(ONE/PQALEN+ONE/PQBLEN)/QLEN/TWOPI
      END IF


C IF LPONEL THEN SOLUTION IS COMPLETE

      IF (LPONEL) GOTO 999


C LOOP THROUGH THE QUADRATURE POINTS AND ACCUMULATE INTEGRALS

        DO 120 IQ=1,NQ
C A: Set Q(1..2) the point on the element.
          Q(1)=QA(1)+AQ(IQ)*QBMA(1)
          Q(2)=QA(2)+AQ(IQ)*QBMA(2)
C B: Set NORMQ(1..2) (already set).
C C: Compute DNPNQ [dot product of VECP and NORMQ] (already set).
C D: Compute RR(1..2) = P(1..2)-Q(1..2)
          CALL VEC2(Q,P,RR)
C E: Compute R [modulus of RR(1..2)], SR [R-squared].
          SR=SSIZE2(RR)
          R=SQRT(SR)
C F: Compute CR [R-cubed] (unnecessary).
C G: Compute RNQ [derivative of R with respect to NORMQ].
          IF (LMORN) RNQ=-DOT2(RR,NORMQ)/R
C H: Compute RNP [derivative of R with respect to VECP].
          IF (LMTORN) RNP=DOT2(RR,VECP)/R
C I: Compute RNPRNQ [RNP*RNQ].
          IF (LN) RNPRNQ=RNP*RNQ
C J: Compute RNPNQ [second derivative of R with respect to
C    VECP and NORMQ].
          IF (LN) RNPNQ=-(DNPNQ+RNPRNQ)/R


          IF (LL) THEN
C N: Compute G, the Green's function.
          G=-OO2PI*LOG(R)
C O: Multiply Lk kernel by weight and add to sum SUML.
          SUML=SUML+WQ(IQ)*DBLE(G)
          END IF
C P: Compute GR, the derivative of the G with respect to R.
          IF (LMSORN) GR=-OO2PI/R
C Q: Compute WGR, the weight multiplied by GR.
          IF (LMORMT) WGR=WQ(IQ)*DBLE(GR)
C R: Compute Mk kernel multiplied by weight and add to sum SUMM.
          IF (LM) SUMM=SUMM+DBLE(WGR)*RNQ
C S: Compute Mkt kernel multiplied by weight and add to sum SUMMT.
          IF (LMT) SUMMT=SUMMT+DBLE(WGR)*RNP
          IF (LN) THEN
C T: Compute GRR, the second derivative of G with respect to R.
            GRR=OO2PI/SR
C U: Multiply Nk kernel by weight and add to sum SUMN
            SUMN=SUMN+
     *         WQ(IQ)*(DBLE(GR)*RNPNQ+DBLE(GRR)*RNPRNQ)
            END IF

120       CONTINUE

999     CONTINUE

        IF (LL) DISL=QLEN*SUML
        IF (LM) DISM=QLEN*SUMM
        IF (LMT) DISMT=QLEN*SUMMT
        IF (LN) DISN=QLEN*SUMN
        IF (LPONEL) THEN
          LM=LMSV
          LMT=LMTSV
        END IF
998     CONTINUE
        LFAIL=LERROR
        END


C======================================================================